sapply(c(
         "rjson", 
         "data.table", 
         "dplyr", 
         "ggplot2", 
         "stringr", 
         "purrr", 
         "foreach", 
         "doParallel", 
         "patchwork", 
         "lme4", 
         "lmerTest",
         "testit",
         "latex2exp",
         "ggpubr"
         ), 
       require, character=TRUE)
## Loading required package: rjson
## Loading required package: data.table
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: stringr
## Loading required package: purrr
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
## 
##     transpose
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: doParallel
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: patchwork
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: lmerTest
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
## Loading required package: testit
## Loading required package: latex2exp
## Loading required package: ggpubr
##      rjson data.table      dplyr    ggplot2    stringr      purrr    foreach 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
## doParallel  patchwork       lme4   lmerTest     testit  latex2exp     ggpubr 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE
sf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()
##         ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value   ?                               ?                                 
## visible FALSE                           FALSE                             
##         ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value   ?                           ?                                          
## visible FALSE                       FALSE                                      
##         ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value   ?                                                  
## visible FALSE                                              
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value   ?                                                                
## visible FALSE                                                            
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/Models/RunMQLearner.R
## value   ?                                      
## visible FALSE                                  
##         ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value   ?                                                            
## visible FALSE                                                        
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value   ?                                                          
## visible FALSE                                                      
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value   ?                                                                       
## visible FALSE                                                                   
##         ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value   ?                                                              
## visible FALSE                                                          
##         ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value   ?                                                      
## visible FALSE                                                  
##         ./Functions/Model/ModelUtils/aModelUtils.R
## value   ?                                         
## visible FALSE                                     
##         ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value   ?                                                             
## visible FALSE                                                         
##         ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/PrepDfForModel.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/ModelUtils/RecodeDfs.R
## value   ?                                       
## visible FALSE                                   
##         ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value   ?                                                        
## visible FALSE                                                    
##         ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value   ?                                          
## visible FALSE                                      
##         ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value   ?                                        
## visible FALSE                                    
##         ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value   ?                                             
## visible FALSE                                         
##         ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value   ?                                                 
## visible FALSE                                             
##         ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value   ?                                               
## visible FALSE                                           
##         ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value   ?                                              
## visible FALSE                                          
##         ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value   ?                                                   
## visible FALSE                                               
##         ./Functions/Model/TrialWiseComps/DecayQVals.R
## value   ?                                            
## visible FALSE                                        
##         ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value   ?                                                
## visible FALSE
DefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))

Load data from studies 1 and 2

s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv")
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv")
s1_pst <- read.csv("../data/cleaned_files/s1_PST_deident.csv") %>% rename(ID=deident_ID)
s2_pst <- read.csv("../data/cleaned_files/s2_PST_deident.csv") %>% rename(ID=deident_ID)

Harmonize variables and create some separate vars

# Study 2 harmonize  
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive" 
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2

s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))

s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))

Define paths and functions

# MLE results  
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results 
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations  
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function 
rm <- function(path, model_str) read.csv(paste0(path, model_str))

Plotting delay

s1_train_delay_summs <- s1_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.
s1_tr_delay_summs <- 
  Rmisc::summarySEwithin(s1_train_delay_summs,
                        measurevar = "m",
                        withinvars = c("condition", "delay_trunc"), 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition, delay_trunc
s2_train_delay_summs <- s2_train %>% filter(!is.na(delay_trunc)) %>% group_by(condition, delay_trunc, ID) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'condition', 'delay_trunc'. You can
## override using the `.groups` argument.
s2_tr_delay_summs <- 
  Rmisc::summarySEwithin(s2_train_delay_summs,
                        measurevar = "m",
                        withinvars = c("condition", "delay_trunc"), 
                        idvar = "ID")
## Automatically converting the following non-factors to factors: condition, delay_trunc
a <- ggplot(s1_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) + 
  geom_line(size=2, position = position_dodge(width = .2)) +
  geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") + 
  geom_hline(yintercept = .5, size=2, color="black") +
  geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
                position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
  geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
  ga + ap + ft + tol + ylim(.48, .82) + tp +
  scale_color_manual(values=c("purple", "orange"), 
                     labels=c("cognitive", "overt"))  + ylab("Proportion correct") + xlab("Delay") + 
  ggtitle("Experiment 1")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
b <- ggplot(s2_tr_delay_summs, aes(x=delay_trunc, y=m, group=condition, color=condition)) + 
  geom_line(size=2, position = position_dodge(width = .2)) +
  geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3, color="gray57") + 
  geom_hline(yintercept = .5, size=2, color="black") +
  geom_errorbar(aes(x=delay_trunc, y=m, ymin=m-se, ymax=m+se, group=condition),
                position = position_dodge(width = .2), color="black", size=1.5, width=.15) +
  geom_point(size=6, position = position_dodge(width = .2), pch=21, fill="gray57", alpha=.95) +
  ga + ap + ft + lp + ylim(.48, .82) + tp +
  scale_color_manual(values=c("purple", "orange"), 
                     labels=c("cognitive", "overt"))  + ylab("") + xlab("Delay") + 
  ggtitle("Experiment 2")
delay_comb_plot <- a+b
delay_comb_plot

#ggsave("../paper/figs/fig_parts/delay_plot.png", delay_comb_plot, height=7, width=10, dpi=700)

Supplemental material model validation (although learning curves exp 2 plot is in key_results)

6/5/23 — reran sims to mkae sure completely updated/bug fixed version of par results used

s1_train_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_train_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")

s1_test_sim_m1_eb <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")

s2_test_sim_m1_eb <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")

Capture of delay effects

out_a <- SimplePlotSimTrainingDelay(emp_df=s1_train, sim_df=s1_train_sim_m1_eb)
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## Warning: Duplicated aesthetics after name standardisation: colour
out_b <- SimplePlotSimTrainingDelay(emp_df=s2_train, sim_df=s2_train_sim_m1_eb) 
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'iter'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.
## Warning: Duplicated aesthetics after name standardisation: colour
out_a
## Warning: Duplicated aesthetics after name standardisation: colour

out_b
## Warning: Duplicated aesthetics after name standardisation: colour

# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_1.png", out_a, height=8, width=12, dpi=700)
# ggsave("../paper/figs/fig_parts/simvsemp_delay_plot_2.png", out_b, height=8, width=12, dpi=700)

Worst-option Proportions in Test Phase and comparison to same model without CK

worst_plots_a_ck <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")
w_a_ck <- worst_plots_a_ck[[1]]+worst_plots_a_ck[[2]]
w_a_ck 

worst_plots_b_ck <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb, m_string="Simulations \nwith Choice Kernel")
w_b_ck <- worst_plots_b_ck[[1]]+worst_plots_b_ck[[2]]
w_b_ck
## Warning: Removed 1 rows containing missing values (`geom_point()`).

# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_1.png", w_a_ck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/CK_simvsemp_worst_plot_2.png", w_b_ck, height=6, width=15, dpi=700)

Comparative worst plots

s1_test_sim_m1_eb_no_ck <- 
  rm(sp, 
     "SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK16278.csv")

s2_test_sim_m1_eb_no_ck <- 
  rm(sp, 
    "SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPriorNoCK31508.csv")
worst_plots_nck_a <- PlotAllWorstTestOptionsSimVsEmp(s1_sit, s1_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")
w_a_nck <- worst_plots_nck_a[[1]]+worst_plots_nck_a[[2]]
w_a_nck

worst_plots_nck_b <- PlotAllWorstTestOptionsSimVsEmp(s2_sit, s2_test_sim_m1_eb_no_ck, m_string="Simulations \nNo Choice Kernel")

Captures the mis-spec

w_b_nck <- worst_plots_nck_b[[1]]+worst_plots_nck_b[[2]]
w_b_nck

# ggsave("../paper/figs/fig_parts/NOCK_simvsemp_worst_plot_1.png", w_a_nck, height=6, width=15, dpi=700)
# ggsave("../paper/figs/fig_parts/NOCKsimvsemp_worst_plot_2.png", w_b_nck, height=6, width=15, dpi=700)

Parameter recovery

m1_par_recov <- 
  read.csv("../model_res/par_recov_clean/pr_opts_mle/best_pr/par_recov_BEST_EMPIRICAL_BAYES_mv_gaussstudy_1_train_SIT_RunMQLearnerDiffDecayToPessPrior_76217.csv") # Confirmed created 3/28 after the 3/23 bug fix on decay  
eb_recov_pars_m1 <- m1_par_recov[grep("EB_recovered", names(m1_par_recov))]
simmed_pars_m1 <- m1_par_recov[grep("simmed", names(m1_par_recov))]

Phi difference

this_eb_pr_df <- 
  data.frame("simulated"=simmed_pars_m1$cog_phi_simmed-simmed_pars_m1$overt_phi_simmed,
             "EB_recovered"=eb_recov_pars_m1$cog_phi_EB_recovered-eb_recov_pars_m1$overt_phi_EB_recovered)

phi_diff <- ggplot(this_eb_pr_df, aes(x=simulated, y=EB_recovered)) +
      geom_line(aes(x=simulated, y=simulated), size=4) +
      geom_point(size=8, pch=21, fill="gray57", alpha=.8) + 
      stat_cor(method="spearman", size=6, label.y=.9) +
      ggtitle(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
      ga + ap + tp  + 
      ylab("Recovered") + xlab("Simulated")
phi_diff

#ggsave("../paper/figs/fig_parts/pr/pr_phi_diff.png", phi_diff, height=4, width=8, dpi=700)

Percent correctly recovering the sign difference

this_eb_pr_df$sim_sign <- if_else(this_eb_pr_df$simulated >= 0, 1, 0)
this_eb_pr_df$eb_sign <- if_else(this_eb_pr_df$EB_recovered >= 0, 1, 0)
# this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign
# (this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1
ta <- table((this_eb_pr_df$sim_sign==this_eb_pr_df$eb_sign)*1)
ta[2]/sum(ta)
##         1 
## 0.7142857
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$epsilon_simmed, eb_recov_pars_m1$epsilon_EB_recovered, 
                     "$\\epsilon", stat_pos = .32)
plot

# ggsave("../paper/figs/fig_parts/pr/eps.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$cog_phi_simmed, eb_recov_pars_m1$cog_phi_EB_recovered, 
                     "$\\phi^{Cog}")
plot

# ggsave("../paper/figs/fig_parts/pr/phi_cog.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$overt_phi_simmed, eb_recov_pars_m1$overt_phi_EB_recovered, 
                     "$\\phi^{Overt}")
plot

# ggsave("../paper/figs/fig_parts/pr/phi_overt.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$choice_LR_simmed, eb_recov_pars_m1$choice_LR_EB_recovered, 
                     "$\\alpha_{CK}")
plot

# ggsave("../paper/figs/fig_parts/pr/ck.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$explor_scalar_simmed, eb_recov_pars_m1$explor_scalar_EB_recovered, 
                     "ES")
plot

# ggsave("../paper/figs/fig_parts/pr/es.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$q_LR_simmed, eb_recov_pars_m1$q_LR_EB_recovered, 
                     "$\\alpha_{Q}")
plot

# ggsave("../paper/figs/fig_parts/pr/qlr.png",
#        plot, height=4, width=8, dpi=700)
plot <- NULL
plot <- CreatePRPlot(simmed_pars_m1$beta_simmed, eb_recov_pars_m1$beta_EB_recovered, 
                     "$\\beta", stat_pos = 25)
plot

Simulation with a high difference in phi

Load model fits

m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")

m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")

m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>% 
  group_by(ID) %>% slice(which.min(nll))

Find a participant with relatively high diff in decay

m1_study1_eb$cog_phi[25]
## [1] 0.6817201
m1_study1_eb$overt_phi[25]
## [1] 0.1300788
m1_study1_eb[25, ]
## # A tibble: 1 × 12
## # Groups:   ID [1]
##       X epsilon  q_LR cog_phi overt_phi choice_LR explor_scalar  beta
##   <int>   <dbl> <dbl>   <dbl>     <dbl>     <dbl>         <dbl> <dbl>
## 1    25  0.0365 0.732   0.682     0.130    0.0425         0.220  4.97
## # … with 4 more variables: convergence <int>, nll <dbl>, AIC <dbl>, ID <int>

This did correspond to a big difference in performance empirically

s1_train %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.
## # A tibble: 8 × 3
## # Groups:   valence_and_probability [4]
##   valence_and_probability condition     m
##   <chr>                   <chr>     <dbl>
## 1 40-10_punishment        cognitive 0.6  
## 2 40-10_punishment        overt     0.95 
## 3 40-10_reward            cognitive 0.9  
## 4 40-10_reward            overt     1    
## 5 90-10_punishment        cognitive 0.525
## 6 90-10_punishment        overt     0.775
## 7 90-10_reward            cognitive 0.05 
## 8 90-10_reward            overt     1
s1_sit %>% filter(ID==25) %>% group_by(valence_and_probability, condition) %>% summarize(m=mean(correct))
## `summarise()` has grouped output by 'valence_and_probability'. You can override
## using the `.groups` argument.
## # A tibble: 8 × 3
## # Groups:   valence_and_probability [4]
##   valence_and_probability condition     m
##   <chr>                   <chr>     <dbl>
## 1 40-10_punishment        cognitive  0.75
## 2 40-10_punishment        overt      1   
## 3 40-10_reward            cognitive  1   
## 4 40-10_reward            overt      1   
## 5 90-10_punishment        cognitive  1   
## 6 90-10_punishment        overt      1   
## 7 90-10_reward            cognitive  0   
## 8 90-10_reward            overt      1

Summarize the correct and incorrect q-values

corr_q_vals_summ <- 
  s1_train_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_correct_Q_vals))
## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.
corr_q_vals_summ$valence <- factor(corr_q_vals_summ$valence, levels=c("reward", "punishment"))
corr_q_vals_summ[corr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"
incorr_q_vals_summ <- 
  s1_train_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_incorrect_Q_vals))
## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.
incorr_q_vals_summ$valence <- factor(incorr_q_vals_summ$valence, levels=c("reward", "punishment"))

incorr_q_vals_summ[incorr_q_vals_summ$cond=="cog", "cond"] <- "cognitive"

Summarize the simulated performance during learning

sim_perf_summ <- 
  s1_train_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence, trial_within_condition) %>% summarize(m=mean(sim_corrs))
## `summarise()` has grouped output by 'cond', 'probability', 'valence'. You can
## override using the `.groups` argument.
sim_perf_summ$valence <- factor(sim_perf_summ$valence, levels=c("reward", "punishment"))

sim_perf_summ[sim_perf_summ$cond=="cog", "cond"] <- "cognitive"

Take the difference in q-values

assert(all(corr_q_vals_summ$trial_within_condition==incorr_q_vals_summ$trial_within_condition))
assert(all(corr_q_vals_summ$valence==incorr_q_vals_summ$valence))
assert(all(corr_q_vals_summ$probability==incorr_q_vals_summ$probability))

corr_q_vals_summ$diff <- corr_q_vals_summ$m-incorr_q_vals_summ$m
a <- 
  ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=m, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.9, color="blue", size=1.5, linetype="longdash") +
  geom_hline(yintercept=-.1, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Q-value") +
  ggtitle("Correct") +
  xlab("Stim iteration") + ylim(-1, 1)
b <- 
  ggplot(incorr_q_vals_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=m, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.1, color="blue", size=1.5, linetype="longdash") +
  geom_hline(yintercept=-.9, color="red", size=1.5, linetype="longdash") +
  ga + ap + tp + ft + tol +
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ggtitle("Incorrect") +
  ylab("") +
  xlab("Stim iteration") + ylim(-1, 1)
c <- 
   ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=diff, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
  geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Difference") +
  xlab("Stim iteration") + ylim(-1, 1)

# For legend  
# ggplot(corr_q_vals_summ %>% filter(probability=="90-10"),   
#        aes(x=trial_within_condition, y=diff, color=valence)) +
#   geom_hline(yintercept=0, color="black", size=2) + 
#   geom_line(size=3) + facet_wrap( ~ cond) + 
#   scale_color_manual(values=c("blue", "red")) +
#   geom_hline(yintercept=.8, color="blue", size=1.5, linetype="longdash") +
#   geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
#   ga + ap +   ft + 
#   theme(legend.text = element_text(size = 30),
#                legend.title = element_blank(),
#                legend.key.size = unit(2.5, 'lines')) +
#   theme(plot.title = element_text(size = 30, hjust = .5)) +
#   ylab("") +
#   ggtitle("Difference") +
#   xlab("Stim iteration") + ylim(-1, 1)
d <- ggplot(sim_perf_summ %>% filter(probability=="90-10"),   
       aes(x=trial_within_condition, y=m, color=valence)) +
  geom_hline(yintercept=0.5, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ cond) + 
  scale_color_manual(values=c("blue", "red")) +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Proportion correct") +
  ggtitle("Learning") +
  xlab("Stim iteration") 
sim_test_perf_summ <- 
  s1_test_sim_m1_eb %>% filter(ID==25) %>% 
  group_by(cond, probability, valence) %>% summarize(m=mean(sit_sim_corrs))
## `summarise()` has grouped output by 'cond', 'probability'. You can override
## using the `.groups` argument.
sim_test_perf_summ$valence <- factor(sim_test_perf_summ$valence, levels=c("reward", "punishment"))

sim_test_perf_summ[sim_test_perf_summ$cond=="cog", "cond"] <- "cognitive"
#sim_test_perf_summ %>% filter(probability=="90-10")
e <- ggplot(sim_test_perf_summ %>% filter(probability=="90-10"),   
       aes(x=valence, y=m, color=valence)) + 
  geom_hline(yintercept=0.5, color="black", size=2) +
  geom_bar(stat="identity", fill="white", size=3) + 
  facet_wrap(~cond) +
  
  scale_color_manual(values=c("blue", "red")) +
  ga + ap + tol + ft +
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Proportion correct") +
  ggtitle("Test") +
  
  xlab("") + theme(axis.ticks.x=element_blank(), axis.text.x=element_blank())#+ ylim(.48, 1) +
q_vals_comb <- a + b + c
q_vals_comb

sim_perf_comb <- d + e
sim_perf_comb

#ggsave("../paper/figs/fig_parts/q_vals_plot.png", q_vals_comb, height=5, width=16, dpi=700)
#ggsave("../paper/figs/fig_parts/sim_perf_comb.png", sim_perf_comb, height=6, width=14, dpi=700)

Supplementary Figure 7

Plotting Q-values by valence for m3 vs. m4 (m6 vs. m7) to understand how the introduction of pessimistic priors affects performance

# Pessimistic  
m3_sims_s1 <-  
  rm(sp,  "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayToPessPrior82512.csv")
# Naive  
m4_sims_s1 <-  
  rm(sp, "SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDecayTo0Inits36394.csv")

M6: Q nothing varied in paper

pess_priors_qv_summs <- m3_sims_s1 %>% filter(probability=="90-10") %>% 
  group_by(valence, trial_within_condition) %>% 
  summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>% 
  mutate(diff=m_corr-m_incorr)
## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.
pess_priors_qv_summs$valence <- factor(pess_priors_qv_summs$valence, levels=c("reward", "punishment"))

M7: Q nothing varied, decay to 0

naive_priors_qv_summs <- m4_sims_s1 %>% filter(probability=="90-10") %>% 
  group_by(valence, trial_within_condition) %>% 
  summarize(m_corr=mean(sim_correct_Q_vals), m_incorr=mean(sim_incorrect_Q_vals)) %>% 
  mutate(diff=m_corr-m_incorr)
## `summarise()` has grouped output by 'valence'. You can override using the
## `.groups` argument.
naive_priors_qv_summs$valence <- factor(naive_priors_qv_summs$valence, levels=c("reward", "punishment"))

Break down into correct/incorrect plots

ggplot(pess_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_corr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Correct") +
  xlab("Stim iteration") + ylim(-1, 1)

ggplot(naive_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_corr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Correct") +
  xlab("Stim iteration") + ylim(-1, 1)

ggplot(pess_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_incorr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  #geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Incorrect") +
  xlab("Stim iteration") + ylim(-1, 1)

ggplot(naive_priors_qv_summs,   
       aes(x=trial_within_condition, y=m_incorr, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  # geom_hline(yintercept=.81, color="red", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Incorrect") +
  xlab("Stim iteration") + ylim(-1, 1)

pess_plot <- ggplot(pess_priors_qv_summs,   
       aes(x=trial_within_condition, y=diff, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("Q-value Difference: \n Correct vs. Incorrect") +
  ggtitle("Pessimistic Priors") +
  xlab("Stim iteration") + ylim(0, 1)
naive_plot <- ggplot(naive_priors_qv_summs,   
       aes(x=trial_within_condition, y=diff, color=valence)) +
  geom_hline(yintercept=0, color="black", size=2) + 
  geom_line(size=3) + facet_wrap( ~ valence) + 
  scale_color_manual(values=c("blue", "red")) +
  geom_hline(yintercept=.8, color="black", size=1.5, linetype="longdash") +
  ga + ap + tol + ft + 
  theme(plot.title = element_text(size = 30, hjust = .5)) +
  ylab("") +
  ggtitle("Naive Priors (All Zero)") +
  xlab("Stim iteration") + ylim(0, 1)
simple_qv_plots_comb <- pess_plot + naive_plot
simple_qv_plots_comb 

No powerpoint — loaded directly into supplemental

#ggsave("../paper/figs/fig_parts/simple_qv_plots_comb.png", simple_qv_plots_comb, height=7, width=14, dpi=700)

PST/Stimulus valuation phase analyses

Used rew history because Q-values are for specific Q(s,a)s whereas selection in the PST phase is between stimuli ie. Q(s, :)

# Take just the trials within condition   
within_cond_pst_s1 <- s1_pst %>% filter(test_condition == "cogn_cogn" | test_condition == "overt_overt")
within_cond_pst_s1$test_condition <- factor(within_cond_pst_s1$test_condition, levels=c("overt_overt", "cogn_cogn"))
contrasts(within_cond_pst_s1$test_condition) <- c(-.5, .5)
head(within_cond_pst_s1$test_condition)
## [1] overt_overt overt_overt overt_overt overt_overt overt_overt overt_overt
## attr(,"contrasts")
##             [,1]
## overt_overt -0.5
## cogn_cogn    0.5
## Levels: overt_overt cogn_cogn
# Both singular  
# summary(m0_s1 <- glmer(resp_num ~  scale(left_min_right) + ( scale(left_min_right)|ID),
#           data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# 
# summary(m0_s1 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (scale(left_min_right) |ID),
#           data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))

Using model without random intercepts because not singular

summary(m0_s1 <- glmer(resp_num ~  scale(left_min_right) + (0 + scale(left_min_right)|ID),
          data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   6456.5   6476.6  -3225.2   6450.5     5983 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3395 -0.7468  0.1149  0.7468  4.5130 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.538    1.24    
## Number of obs: 5986, groups:  ID, 125
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.11210    0.03073   3.648 0.000264 ***
## scale(left_min_right)  1.41202    0.12084  11.685  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(lft_m_) 0.015
summary(m1_s1 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
          data=within_cond_pst_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s1
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   6431.2   6464.7  -3210.6   6421.2     5981 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.1380 -0.7457  0.1136  0.7448  4.1139 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.555    1.247   
## Number of obs: 5986, groups:  ID, 125
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            0.11079    0.03085   3.592 0.000328 ***
## scale(left_min_right)                  1.42529    0.12150  11.731  < 2e-16 ***
## test_condition1                        0.09218    0.06167   1.495 0.134955    
## scale(left_min_right):test_condition1 -0.36515    0.07029  -5.195 2.05e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sc(__) tst_c1
## scl(lft_m_)  0.015              
## test_cndtn1 -0.041  0.005       
## scl(l__):_1  0.003 -0.035  0.025
# Take just the trials within condition   
within_cond_pst_s2 <- s2_pst %>% filter(test_condition == "cognitive_cognitive" | test_condition == "overt_overt")
within_cond_pst_s2$test_condition <- 
  factor(within_cond_pst_s2$test_condition, levels=c("overt_overt", "cognitive_cognitive"))
contrasts(within_cond_pst_s2$test_condition) <- c(-.5, .5)
head(within_cond_pst_s2$test_condition)
## [1] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## [4] cognitive_cognitive cognitive_cognitive cognitive_cognitive
## attr(,"contrasts")
##                     [,1]
## overt_overt         -0.5
## cognitive_cognitive  0.5
## Levels: overt_overt cognitive_cognitive

Use the comparable model

summary(m0_s2 <- glmer(resp_num ~  scale(left_min_right) + (0 + scale(left_min_right)|ID),
          data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: resp_num ~ scale(left_min_right) + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7161.3   7181.7  -3577.7   7155.3     6616 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.1402 -0.7899  0.1255  0.7803  4.7614 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.507    1.228   
## Number of obs: 6619, groups:  ID, 138
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.11126    0.02914   3.818 0.000134 ***
## scale(left_min_right)  1.35153    0.11345  11.912  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(lft_m_) 0.014
summary(m1_s2 <- glmer(resp_num ~  scale(left_min_right)*test_condition + (0 + scale(left_min_right) |ID),
          data=within_cond_pst_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## resp_num ~ scale(left_min_right) * test_condition + (0 + scale(left_min_right) |  
##     ID)
##    Data: within_cond_pst_s2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   7163.5   7197.5  -3576.8   7153.5     6614 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.0975 -0.7881  0.1258  0.7749  4.9593 
## 
## Random effects:
##  Groups Name                  Variance Std.Dev.
##  ID     scale(left_min_right) 1.509    1.228   
## Number of obs: 6619, groups:  ID, 138
## 
## Fixed effects:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            0.11101    0.02914   3.809  0.00014 ***
## scale(left_min_right)                  1.35180    0.11350  11.910  < 2e-16 ***
## test_condition1                        0.06205    0.05826   1.065  0.28689    
## scale(left_min_right):test_condition1 -0.05321    0.06672  -0.797  0.42518    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sc(__) tst_c1
## scl(lft_m_)  0.014              
## test_cndtn1 -0.008  0.004       
## scl(l__):_1  0.006 -0.002  0.027
p_s1 <- 
  sjPlot::plot_model(m1_s1, type="pred", terms = c("left_min_right [all]", "test_condition"))


p_fin_1 <- p_s1 + ap + tp + ga + 
  xlab("Left minus right reward history") + ylab("") + ggtitle("") + 
  theme(plot.title = element_text(hjust=0)) +
  ylab("Chance of picking left") +
   xlab("Left minus right reward history") +
  ggtitle("Predicted probabilities") + 
  scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) + tol
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
p_s2 <- 
  sjPlot::plot_model(m1_s2, type="pred", terms = c("left_min_right [all]", "test_condition"))


p_fin_s2 <- p_s2 + ap + tp + ga + 
  xlab("Left minus right reward history") + ylab("") + ggtitle("") + 
  theme(plot.title = element_text(hjust=0)) +
  ylab("Chance of picking left") +
   xlab("Left minus right reward history") +
  #ggtitle("Predicted probabilities") + 
  scale_color_manual(values=c("orange", "purple"), labels=c("overt-overt", "cognitive-cognitive")) +
  theme(legend.text = element_text(size = 18),
               legend.title = element_blank(),
               legend.key.size = unit(1.3, 'lines'))
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
p_fin_1 + p_fin_s2

Mixed test phase

unique(s1_pst$test_condition)
## [1] "overt_cogn"  "overt_overt" "cogn_cogn"
mix_s1 <- s1_pst %>% filter(test_condition == "overt_cogn")
# If cognitive is on the left and they chose left then they chose cognitive  
mix_s1[mix_s1$left_training_cond=="cognTraining", "chose_cognitive"] <- 
  if_else(mix_s1[mix_s1$left_training_cond=="cognTraining", "response"]=="left", 1, 0)

# If cognitive is on the right and they chose right then they chose cognitive
mix_s1[mix_s1$right_training_cond=="cognTraining", "chose_cognitive"] <- 
  if_else(mix_s1[mix_s1$right_training_cond=="cognTraining", "response"]=="right", 1, 0)
table(mix_s1$chose_cognitive)[2]/sum(table(mix_s1$chose_cognitive))
##         1 
## 0.4751506
unique(s2_pst$test_condition)
## [1] "overt_cognitive"     "cognitive_cognitive" "overt_overt"
mix_s2 <- s2_pst %>% filter(test_condition == "overt_cognitive")
# If cognitive is on the left and they chose left then they chose cognitive  
mix_s2[mix_s2$left_training_cond=="cognitive", "chose_cognitive"] <- 
  if_else(mix_s2[mix_s2$left_training_cond=="cognitive", "response"]=="left", 1, 0)

# If cognitive is on the right and they chose right then they chose cognitive
mix_s2[mix_s2$right_training_cond=="cognitive", "chose_cognitive"] <- 
  if_else(mix_s2[mix_s2$right_training_cond=="cognitive", "response"]=="right", 1, 0)
table(mix_s2$chose_cognitive)[2]/sum(table(mix_s2$chose_cognitive))
##         1 
## 0.4231989
mixed_summs_s1 <- mix_s1 %>% group_by(choice_type_plotting) %>% 
  summarize(m=mean(chose_cognitive), n())
mixed_summs_s1
## # A tibble: 4 × 3
##   choice_type_plotting     m `n()`
##   <chr>                <dbl> <int>
## 1 "freq-P\nfreq-P"     0.490   997
## 2 "freq-R\nfreq-R"     0.428   997
## 3 "infreq-P\ninfreq-P" 0.553   992
## 4 "infreq-R\ninfreq-R" 0.429   998
mixed_summs_s2 <- mix_s2 %>% group_by(choice_type_plotting) %>% 
  summarize(m=mean(chose_cognitive), n())
mixed_summs_s2
## # A tibble: 4 × 3
##   choice_type_plotting     m `n()`
##   <chr>                <dbl> <int>
## 1 "freq-P\nfreq-P"     0.408  1104
## 2 "freq-R\nfreq-R"     0.393  1104
## 3 "infreq-P\ninfreq-P" 0.447  1102
## 4 "infreq-R\ninfreq-R" 0.445  1104
summary(choice_effect_s1 <- glmer(chose_cognitive ~  1 + ( 1 |ID),
          data=mix_s1, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
##    Data: mix_s1
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   5269.0   5281.6  -2632.5   5265.0     3982 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4732 -0.8937 -0.4973  0.9600  2.3488 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.5492   0.7411  
## Number of obs: 3984, groups:  ID, 125
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.10926    0.07443  -1.468    0.142
summary(choice_effect_s2 <- glmer(chose_cognitive ~  1 + ( 1|ID),
          data=mix_s2, family="binomial", control = glmerControl(optimizer = "bobyqa")))
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: chose_cognitive ~ 1 + (1 | ID)
##    Data: mix_s2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   5829.4   5842.2  -2912.7   5825.4     4412 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5760 -0.7917 -0.5754  1.0417  2.2498 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.3916   0.6258  
## Number of obs: 4414, groups:  ID, 138
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.33769    0.06215  -5.433 5.53e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Questionnaires

s1_qs <- read.csv("../data/cleaned_files/s1_qs_deidentified.csv")
s2_qs <- read.csv("../data/cleaned_files/s2_qs_deidentified.csv")

Exclude pts who answered both catch questions incorrectly

S1 all correct so no exclusions

table(s1_qs$Q6_20) # All correct
## 
## A little bit 
##          123
table(s1_qs$Q7_19) # All correct
## 
## Very false for me 
##               123
table(s2_qs$Q6_20) # All correct
## 
## A little bit 
##          137
table(s2_qs$Q7_19) # One pt incorrect but they answered the other catch item correctly, so retained (given recs in lit that excluding on a single catch item can bias results)
## 
## Very false for me  Very true for me 
##               136                 1
#s2_qs %>% filter(Q7_19 == "Very true for me")

S1 questionnaires

Get and recode PTQ and RRS

s1_ptq <- s1_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
s1_ptq$ID <- s1_qs$ID
s1_ptq_copy <- s1_ptq
s1_ptq_copy$ID <- s1_qs$ID
s1_ptq[s1_ptq=="Never"] <- 0
s1_ptq[s1_ptq=="Rarely"] <- 1
s1_ptq[s1_ptq=="Sometimes"] <- 2 
s1_ptq[s1_ptq=="Often"] <- 3
s1_ptq[s1_ptq=="Almost always"] <- 4
s1_rrs <- 
  s1_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
s1_rrs$ID <- s1_qs$ID
s1_rrs_copy <- s1_rrs
s1_rrs_copy$ID <- s1_qs$ID
s1_rrs[s1_rrs=="Almost never"] <- 1
s1_rrs[s1_rrs=="Sometimes"] <- 2
s1_rrs[s1_rrs=="Often"] <- 3
s1_rrs[s1_rrs=="Almost always"] <- 4

Drop the 1 pt who skipped all items (leaving n=122 completers for study 1)

s1_ptq <- s1_ptq[-66, ]
s1_rrs <- s1_rrs[-66, ]

Convert to numeric

s1_ptq_num <- foreach (i = 1:ncol(s1_ptq)) %do% {
  data.frame(as.numeric(unlist(s1_ptq[i]))) %>% setNames(names(s1_ptq[i]))
}%>% bind_cols()

Just 4 data points (< .1% of data) missing for remaining, so mean impute these

length(which(is.na(s1_ptq_num[, 1:15])))
## [1] 4
length(which(is.na(s1_ptq_num[, 1:15])))/(dim(s1_ptq_num[, 1:15])[1]*dim(s1_ptq_num[, 1:15])[2]) 
## [1] 0.002185792
s1_ptq_final <- foreach (i = 1:nrow(s1_ptq_num)) %do% {
  if (any(is.na(s1_ptq_num[i, ]))) {
    s1_ptq_num[i, is.na(s1_ptq_num[i, ])] <- mean(unlist(s1_ptq_num[i, ]), na.rm=TRUE)
  }
s1_ptq_num[i, ]
}%>% bind_rows()
hist(rowSums(s1_ptq_final[, 1:15]), breaks=25)

s1_ptq_final$ptq_sum <- rowSums(s1_ptq_final[, 1:15])

Spot check IDs lined up properly

# s1_ptq_final %>% filter(ID == 22)
# s1_ptq_copy %>% filter(ID == 22)
# s1_qs %>% filter(ID == 22)
# s1_ptq_final %>% filter(ID == 104)
# s1_ptq_copy %>% filter(ID == 104)
# s1_qs %>% filter(ID == 104)

Convert to numeric

s1_rrs_num <- foreach (i = 1:ncol(s1_rrs)) %do% {
  data.frame(as.numeric(unlist(s1_rrs[i]))) %>% setNames(names(s1_rrs[i]))
}%>% bind_cols()
length(which(is.na(s1_rrs_num[, 1:10])))
## [1] 3
length(which(is.na(s1_rrs_num[, 1:10])))/(dim(s1_rrs_num[, 1:10])[1]*dim(s1_rrs_num[, 1:10])[2]) 
## [1] 0.002459016
s1_rrs_final <- foreach (i = 1:nrow(s1_rrs_num)) %do% {
  if (any(is.na(s1_rrs_num[i, ]))) {
    s1_rrs_num[i, is.na(s1_rrs_num[i, ])] <- mean(unlist(s1_rrs_num[i, ]), na.rm=TRUE)
  }
s1_rrs_num[i, ]
}%>% bind_rows()

Spot check

# s1_rrs_final %>% filter(ID == 106)
# s1_rrs_copy %>% filter(ID == 106)
# s1_qs %>% filter(ID == 106)
hist(rowSums(s1_rrs_final[, 1:10]), breaks=25)

s1_rrs_final$rrs_sum <- rowSums(s1_rrs_final[, 1:10])

Sanity check RRS and PTQ are strongly correlated

ComparePars(s1_rrs_final$rrs_sum, s1_ptq_final$ptq_sum, use_identity_line = 0)

Reduced to just matching IDs

m1_s1_model_red <- m1_study1_eb %>% filter(ID %in% s1_ptq_final$ID)
m1_s1_model_red$c_min_o_phi <- m1_s1_model_red$cog_phi - m1_s1_model_red$overt_phi
assert(all(m1_s1_model_red$ID==s1_ptq_final$ID))
assert(all(m1_s1_model_red$ID==s1_rrs_final$ID))

PTQ

ComparePars(m1_s1_model_red$c_min_o_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s1_model_red$cog_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s1_model_red$overt_phi, s1_ptq_final$ptq_sum, use_identity_line = 0)

Same basic pattern for RRS

ComparePars(m1_s1_model_red$c_min_o_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)

# Would expect higher decay to positively correlate w more brooding — instead small and ns neg  
ComparePars(m1_s1_model_red$cog_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)

ComparePars(m1_s1_model_red$overt_phi, s1_rrs_final$rrs_sum, use_identity_line = 0)

s2 questionnaires

Get and recode PTQ and RRS

s2_ptq <- s2_qs %>% select(contains("Q2_")) %>% setNames(paste0("PTQ_", 1:15))
s2_ptq$ID <- s2_qs$ID
s2_ptq_copy <- s2_ptq
s2_ptq_copy$ID <- s2_qs$ID
s2_ptq[s2_ptq=="Never"] <- 0
s2_ptq[s2_ptq=="Rarely"] <- 1
s2_ptq[s2_ptq=="Sometimes"] <- 2 
s2_ptq[s2_ptq=="Often"] <- 3
s2_ptq[s2_ptq=="Almost always"] <- 4
s2_rrs <- 
  s2_qs %>% select(contains("Q8_")) %>% setNames(paste0("RRS_", 1:10))
s2_rrs$ID <- s2_qs$ID
s2_rrs_copy <- s2_rrs
s2_rrs_copy$ID <- s2_qs$ID
s2_rrs[s2_rrs=="Almost never"] <- 1
s2_rrs[s2_rrs=="Sometimes"] <- 2
s2_rrs[s2_rrs=="Often"] <- 3
s2_rrs[s2_rrs=="Almost always"] <- 4

No pts with skips

Convert to numeric

s2_ptq_num <- foreach (i = 1:ncol(s2_ptq)) %do% {
  data.frame(as.numeric(unlist(s2_ptq[i]))) %>% setNames(names(s2_ptq[i]))
}%>% bind_cols()
length(which(is.na(s2_ptq_num[, 1:15])))
## [1] 3
length(which(is.na(s2_ptq_num[, 1:15])))/(dim(s2_ptq_num[, 1:15])[1]*dim(s2_ptq_num[, 1:15])[2]) 
## [1] 0.001459854
s2_ptq_final <- foreach (i = 1:nrow(s2_ptq_num)) %do% {
  if (any(is.na(s2_ptq_num[i, ]))) {
    s2_ptq_num[i, is.na(s2_ptq_num[i, ])] <- mean(unlist(s2_ptq_num[i, ]), na.rm=TRUE)
  }
s2_ptq_num[i, ]
}%>% bind_rows()
hist(rowSums(s2_ptq_final[, 1:15]), breaks=25)

s2_ptq_final$ptq_sum <- rowSums(s2_ptq_final[, 1:15])

Spot check IDs lined up properly

# s2_ptq_final %>% filter(ID == 95)
# s2_ptq_copy %>% filter(ID == 95)
# s2_qs %>% filter(ID == 95)

Convert to numeric

s2_rrs_num <- foreach (i = 1:ncol(s2_rrs)) %do% {
  data.frame(as.numeric(unlist(s2_rrs[i]))) %>% setNames(names(s2_rrs[i]))
}%>% bind_cols()
length(which(is.na(s2_rrs_num[, 1:10])))
## [1] 3
length(which(is.na(s2_rrs_num[, 1:10])))/(dim(s2_rrs_num[, 1:10])[1]*dim(s2_rrs_num[, 1:10])[2]) 
## [1] 0.002189781
s2_rrs_final <- foreach (i = 1:nrow(s2_rrs_num)) %do% {
  if (any(is.na(s2_rrs_num[i, ]))) {
    s2_rrs_num[i, is.na(s2_rrs_num[i, ])] <- mean(unlist(s2_rrs_num[i, ]), na.rm=TRUE)
  }
s2_rrs_num[i, ]
}%>% bind_rows()

Spot check

# s2_rrs_final %>% filter(ID == 128)
# s2_rrs_copy %>% filter(ID == 128)
# s2_qs %>% filter(ID == 128)
hist(rowSums(s2_rrs_final[, 1:10]), breaks=25)

s2_rrs_final$rrs_sum <- rowSums(s2_rrs_final[, 1:10])

Sanity check RRS and PTQ are strongly correlated

ComparePars(s2_rrs_final$rrs_sum, s2_ptq_final$ptq_sum, use_identity_line = 0)

Reduced to just matching IDs

m1_s2_model_red <- m1_study2_eb %>% filter(ID %in% s2_ptq_final$ID)
m1_s2_model_red$c_min_o_phi <- m1_s2_model_red$cog_phi - m1_s2_model_red$overt_phi
assert(all(m1_s2_model_red$ID==s2_ptq_final$ID))
assert(all(m1_s2_model_red$ID==s2_rrs_final$ID))

PTQ

ComparePars(m1_s2_model_red$c_min_o_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s2_model_red$cog_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)

ComparePars(m1_s2_model_red$overt_phi, s2_ptq_final$ptq_sum, use_identity_line = 0)

Same basic pattern for RRS

ComparePars(m1_s2_model_red$c_min_o_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)

# Would expect higher decay to positively correlate w more brooding — instead small and ns neg  
ComparePars(m1_s2_model_red$cog_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)

ComparePars(m1_s2_model_red$overt_phi, s2_rrs_final$rrs_sum, use_identity_line = 0)